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Abstract 

This  study  presents  a  IEEE  standard  Very  High  Speed  Integrated  Circuit  Hardware  Description  Language- Analog  and  Mixed-Signal  Extension 
(VHDL-AMS)  modelling  of  a  complex  multi-domain  energy  conversion  system:  a  fuel  cell  stack.  A  comparative  study  between  the  different 
modelling  approaches  (bond  graphs,  electrical  equivalent  circuits)  is  given  to  show  the  great  advantages  of  the  VHDL-AMS  language  in  the  design 
process  of  fuel  cell  systems.  The  modelling  approach  allows  the  design  team  to  split  the  work  into  several  parts  (concurrent  engineering)  and 
validate  each  part  independently.  The  fuel  cell  stack  model  fits  the  experimental  results.  It  is  able  to  predict  the  voltage  and  the  power  of  the  fuel 
cell  with  a  good  accuracy  taking  into  account  the  water  content  of  the  membrane.  This  last  point  is  really  important  to  design  the  air  supply  system 
(compressor  and  humidifier)  and  its  associated  control. 
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1.  Introduction 

Contemporary  systems  have  become  so  complex  that  mod¬ 
elling  became  a  part  of  the  design  process/cycle  in  order  to 
validate  functionality  and  performance  of  the  whole  system. 
For  this  purpose,  system  designers  need  a  high-level  perspec¬ 
tive  of  the  system.  Fuel  cells,  or  even  more,  fuel  cell  systems 
(FCS)  are  a  good  example  of  complexity,  they  include  multi¬ 
domains  continuous-time  functionalities  like  electromechanical 
(e.g.,  electrical  drives),  electrical  (e.g.,  power  converter,  electri¬ 
cal  loads),  electrochemical,  fluidic  (e.g.,  channels,  compressor) 
and  thermal  parts  [1].  They  also  include  discrete-time  function¬ 
alities  for  the  control  of  such  a  complex  system.  A  fuel  cell 
system  covers  a  range  of  modelling  in  different  energy  domains 
and  spans  levels  of  abstraction  from  low-level  devices  that  make 
the  components  to  the  top-level  functional  unit.  If  we  encompass 
such  range  of  views  of  digital,  analog  and  mixed-signal  systems, 
the  complexity  we  are  dealing  with  is  really  high  and  it  is  not 
possible  to  comprehend  such  complex  systems  in  their  entirety 
without  a  proper  methodology. 
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Methods  which  deal  with  the  complexity  have  to  be  found  in 
order  to  design,  with  some  degree  of  confidence,  components 
and  systems  that  meet  the  requirements.  Systematic  methodol¬ 
ogy  of  design  using  a  top-down  design  approach  has  to  be  applied 
[2] .  This  methodology  decomposes  the  system  design  in  a  collec¬ 
tion  of  components  that  interact  which  can  be  decomposed  until 
a  level  where  we  have  sufficient  details  [3] .  Each  component  can 
have  different  levels  of  details  (compromise  between  accuracy 
and  performance)  and  abstraction.  Each  subsystem  can  be  tested 
as  part  of  the  whole  system  and  can  be  designed  independently 
of  others  what  facilitates  concurrent  engineering  (e.g.,  one  sub¬ 
system  can  be  designed  in  one  place,  the  other  somewhere  else) 
and  improves  productivity. 

For  this  purpose  Hardware  Description  Languages  (HDLs) 
have  been  developed  first  for  the  electronic  (digital)  domain. 
Later  on,  Mixed-Signal  Hardware  Description  Languages 
(MSHDLs)  were  created  to  be  able  to  build  complex  models 
mixing  analog  and  digital  functionalities.  These  programming 
languages  are  used  to  develop  executable  simulation  models 
of  hardware  systems,  not  only  electrical  systems  but  hetero¬ 
geneous  systems  from  several  energy  domains  also  called 
multidisciplinary  or  mixed-disciplines  models  [4].  The  range 
of  applications  resulting  from  the  capability  of  these  powerful 
tools  is  astounding. 
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Modern  HDLs  support  both  behavioural  and  structural  mech¬ 
anisms  [3].  The  first  one  allows  the  designer  to  express  the 
operation  of  a  subsystem  at  various  levels  of  abstraction:  from 
high  abstract  to  very  detailed  levels.  The  second  one,  allows 
the  designer  to  compose  the  model  of  a  complete  system  from 
reusable  model  components  stored  in  a  library.  Thus,  in  top- 
down  design,  an  entire  system,  can  contain  some  components 
with  a  high  level  of  abstraction  described  behaviourally  with 
little  detail,  another  one  can  add,  for  example,  parasitics  and 
another  one  can  decompose  the  component  into  lower  level  ones. 

MSHDLs  support  both  digital  and  analog  functionalities. 
Digital  functionality  allows  to  model  event-driven  techniques 
and  discrete  time.  Analog  functionality  allows  to  model  sys¬ 
tems  of  differential  and  algebraic  equations  (DAEs)  where  the 
solution  varies  continuously  with  time. 

The  needs  for  a  mixed-signal  simulation  capability  have 
resulted  in  the  development  of  several  proprietary  languages 
like  MAST  by  Analogy  in  1986  and  HLD-A  by  Mentor  Graph¬ 
ics  in  1992.  The  languages  of  these  simulation  environments 
are  completely  independent  and  the  use  and  inter-operability 
are  highly  restricted:  models  developed  for  one  system  are  not 
usable  in  other  environments.  Efforts  to  overcome  these  lim¬ 
itations  have  resulted  in  the  development  and  standardisation 
of  VHDL-AMS  (IEEE  standard  Very  High  Speed  Integrated 
Circuit  Hardware  Description  Language-Analog  and  Mixed- 
Signal  Extensions)  presented  hereafter.  Another  open  language, 
Verilog-AMS  (Open  Verilog  International)  has  to  be  mentioned 
but  it  did  not  receive  standardisation  yet.  Interested  readers  can 
refer  to  [4]  where  a  comparison  between  Verilog-AMS  and 
VHDL-AMS  is  given  in  detail. 

This  article  presents  a  VHDL-AMS  model  of  a  fuel  cell  stack. 
In  the  first  part,  the  multi-domain  modelling  language,  VHDL- 
AMS,  is  presented  to  introduce  the  language  basic  principles: 
interested  readers  will  be  able  to  find  some  good  references  to  go 
forward.  In  a  second  part,  a  VHDL-AMS  fuel  cell  stack  model 
is  presented  after  a  comparison  between  the  different  modelling 
approaches  which  have  been  used  in  a  great  deal  of  papers  (bond 
graphs,  electrical  equivalent  circuits  and  signal  flow).  Finally, 
simulation  and  experimental  results  are  presented:  they  show 
very  good  agreements.  The  last  part  concludes  the  work  and 
gives  the  objectives  and  a  roadmap  to  improve  the  model  and 
build  a  Functional  Virtual  Prototype  for  the  design  and  study  of 
fuel  cell  systems. 

2.  Analysis  and  modelling 

2.7.  VHDL-AMS  presentation 
2.1.1.  Overview 

The  IEEE  1076.1  language,  informally  known  as  VHDL- 
AMS,  is  a  superset  of  IEEE  Std  1076-1993  (VHDL)  that 
provides  capabilities  for  describing  and  simulating  analog  and 
mixed- signal  systems  with  conservative  and  nonconservative 
semantics  for  the  analog  portion  of  the  system.  The  language 
supports  many  abstraction  levels  in  electrical  and  nonelectrical 
energy  domains.  The  modelled  analog  systems  are  lumped  sys¬ 
tems  that  can  be  described  by  ordinary  DAEs.  The  language 


does  not  specify  any  particular  technique  to  solve  the  equations, 
but  it  rather  defines  the  results  that  must  be  achieved.  The  solu¬ 
tion  of  the  equations  may  include  discontinuities.  Interaction 
between  the  digital  part  of  a  model  and  its  analog  part  is  sup¬ 
ported  in  a  flexible  and  efficient  manner.  Finally,  support  for 
frequency  domain  small-signal  and  noise  simulation  is  provided 

[5]. 

The  VHDL-AMS  is  designed  to  fill  a  number  of  needs  in  the 
process  design  [6] : 

•  it  allows  description  of  the  structure  of  a  system,  that  is,  how 
it  is  decomposed  into  subsystems  from  different  disciplines 
and  how  those  subsystems  are  interconnected, 

•  it  allows  the  specification  of  the  function  of  a  system  using 
familiar  programming  language  and  equation  forms, 

•  it  allows  the  design  of  a  system  to  be  simulated  before  being 
manufactured:  designers  can  compare  alternatives  and  test 
for  correctness  without  the  delay  and  expense  of  hardware 
prototyping, 

•  it  allows  the  detailed  structure  of  a  design  to  be  synthesised 
from  a  more  abstract  specification,  allowing  designers  to  con¬ 
centrate  on  more  strategic  design  decisions  and  reducing  time 
to  market. 

Moreover,  the  language  is  suited  to  express  plant  models  and 
control  algorithms  needed  for  car  manufacturers  and  suppliers 
[7].  In  other  words,  it  is  able  to  communicate  manufacturer’s 
requests  and  supplier’s  solutions  and  it  is  suitable  for  model 
exchange  within  companies  or  between  researchers  and  compa¬ 
nies. 

Several  commercial  simulators  are  available  today  like  Sys¬ 
tem  Vision  (Mentor  Graphics),  Smash  (Dolphin  Integration), 
Simplorer  (Ansoft)  and  SaberHDL  (Synopsys). 

2.1.2.  Full  working  examples 

In  order  to  illustrate  what  is  a  VHDL-AMS  language,  a  model 
of  a  resistor  will  be  introduced.  Using  VHDL-AMS  terminology, 
the  module,  or  design  entity  is  called  resistor  and  the  input 
and  output  are  ports.  Fig.  1  shows  a  VHDL-AMS  description 
of  the  interface  entity.  This  is  an  example  of  an  entity  decla¬ 
ration.  The  two  first  lines  of  the  listing  indicate  that  the  entity 
involves  of  analog  electrical-energy  systems  package.  Then,  the 
list  of  properties  and  parameters  are  defined.  The  resistor  has  one 
generic  parameter,  that  is  its  resistance  value.  A  default  value  of 
each  generic  parameter  can  be  defined  if  it  is  required.  The  text 
after  are  the  comments  and  are  not  interpreted  by  the  com¬ 
piler.  The  port  names  1 1  and  t2  (line  1 1)  are  analog  terminals 
port  of  the  electrical  nature,  representing  circuit  nodes.  A  graph¬ 
ical  representation  of  the  resistor  entity  is  given  in  Fig.  2.  The 
electrical  nature  specifies  that  the  terminal  has  voltage  and  cur¬ 
rent  properties  associated  with  it.  This  permits  to  connect  the 
resistor  to  other  terminals  which  are  also  of  the  electrical  nature 
type. 

The  next  step  in  building  the  model  is  to  define  the 
behavioural  model  of  the  resistor  as  shown  in  Fig.  3.  Each  ter¬ 
minal  has  to  be  associated  with  across  and  through  branch 
quantities.  This  is  done  in  line  2.  For  electrical  nature,  the  pre- 
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Fig.  1.  VHDL-AMS  resistor  entity  (interface)  declaration. 

(through  quantity) 

ti » -  >  — \A/\A  »t2 

v 

(across  quantity) 

Fig.  2.  Graphical  representation  of  the  VHDL-AMS  resistor  entity. 


library  ieee; 

use  ieee. electrical_systems. all; 

entity  resistor  is 

—  generic  model  parameters 
generic ( 

quantity  R  :  resistance  :=  1.0; 
); 

—  ports  declaration 
port  ( 

terminal  tl,t2  :  electrical 
); 

end  entity  resisitor; 


defined  across  type  is  voltage  and  the  predefined  through  type 
is  current. 

Finally,  a  set  of  differential  and  algebraic  equations  (DAEs) 
can  be  written  to  describe  the  relations  between  the  quantities. 
This  is  done  with  a  simultaneous  statement  (“==”)  which  has 
to  be  verified  at  each  step  of  the  simulation  time.  For  the  resistor, 
the  simultaneous  statement  is  straight  forward:  it  is  the  ohmic 
law  written  in  line  4. 


1  architecture  behav  of  resistor  is 

2  quantity  v  across  i  through  tl  to  t2; 

3  begin 

4  v  ==  i*R; 

5  end  architecture  behav; 
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library  ieee; 

use  ieee. mechanical_systems. all; 
use  ieee. electrical_systems. all; 
use  ieee. math_real. all; 

entity  dcmp  is 
generic ( 

—  armature  resistance 
R  :  resistance  :=  1.0; 

—  armature  inductance 

L  :  inductance  :=  0.01; 

—  motor  torque  constant 
K  :  real  :=  1.0; 

—  armature’s  moment  of  intertia 
J  :  moment _ inertia  :=  0.075; 
-  armature’s  damping 

B  :  damping  :=  0.0); 
port  ( 

—  electrical  terminals 
terminal  tl,  t2  :  electrical; 

—  mechanical  terminal 

—  rotational  velocity  nature 
terminal  m  :  rotational_v; 

—  motor  rotational  speed  in  rpm 
quantity  N  :  out  real; 

); 

Fig.  4.  VHDL-AMS  permanent  magnet  DC  machine. 


Fig.  3.  VHDL-AMS  resistor  architecture  (behaviour)  declaration. 
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end  entity  dcmp; 

architecture  behav  of  dcmp  is 

constant  om2n  :  real:=60.0/(2.0*math_pi) ; 
quantity  v  across  i  through  tl  to  t2; 
quantity  Omega_m  across  T_m  through  m; 
begin 

—  elec.  eq.  :  V  =  Ri  +  +  K Qm 

v  ==  R*i+  L+i’dot+K+Omega.m; 

—  mech.  eq  :  Tm  =  Ki  —  —  B 

T_m  ==  K*i-J*Omega_m,dot“B*Omega_m; 

N  ==  om2n*0mega_m; 

end  architecture  behav; 


It  has  to  be  noticed  that  the  model  does  not  assume  which 
device  is  connected  to  the  machine:  the  causes  why  the  machine 
are  rotating  is  not  known  a  priori.  Consequently,  the  model 
of  the  machine  works  in  the  four  quadrants  depending  on  if 
the  voltage/current  is  imposed  by  a  power  converter  or  if  the 
speed/torque  is  imposed  by  a  mechanical  load. 

2.2.  Fuel  cell  model 

2.2.1.  Fuel  modelling  approach  review 

Several  modelling  approaches  have  been  used  to  model  fuel 
cells:  bond  graphs,  electrical  equivalent  circuits  and  signal  flow. 
One  of  the  first  fuel  cell  multi-domain  approach  has  been 
presented  by  Bernardi  and  Verbrugge  [8]  where  the  authors 
develop  a  one-dimensional  steady  state  and  isothermal  fuel  cell 
model.  The  model  describes  and  predicts  water  transport,  reac¬ 
tant  species  transport  and  ohmic  and  activation  overpotential. 
Multi-domain  approach  involves  derivation  of  a  set  of  equations 
for  each  region  of  the  fuel  cell,  that  is  gas  diffusion  regions 
(anode  and  cathode),  gas  flow  channels,  membrane  and  catalyst 
layer  (anode  and  cathode):  equations  are  solved  separately  and 
simultaneously  [9] .  A  review  and  comparison  of  approaches  for 
modelling  PEM  fuel  cells  have  been  done  in  [9,10]:  interested 
readers  can  refer  to  these  papers  for  more  details.  The  aim  of 
this  part  is  to  compare  the  three  different  approaches  which  are 
bond  graphs,  electrical  equivalent  circuit  and  signal  flow. 


Fig.  4.  ( Continued ). 

More  complex  models  can  be  written  and  not  only  in  the 
electrical  domain.  Fig.  4  shows  a  model  of  a  permanent  mag¬ 
net  DC  machine  which  is  an  electro-mechanical  device  and  has 
consequently  both  electrical  and  mechanical  nature  terminals. 
It  has  two  electrical  terminals  (line  20)  and  one  mechanical  ter¬ 
minal  (shaft)  with  a  rotational- velocity  nature  (line  23).  At  the 
line  25,  a  type  of  port  quantity  is  defined,  the  rotational  speed 
of  the  motor  defined  as  a  out  quantity.  Quantities  can  either  be 
declared  as  in  or  out .  A  graphical  representation  of  the  resistor 
entity  is  given  in  Fig.  5. 

Simultaneous  statements  (lines  34-39)  contain  differential 
equations  which  are  the  well-known  electrical  and  mechanical 
equations  of  a  permanent  DC  machine.  The  time  derivative  of  x 
is  noted  x'  dot. 


Fig.  5.  Graphical  representation  of  the  VHDL-AMS  DC  machine  entity. 


2. 2. 1.1.  Bond  graph  approach.  Bond  graphs  are  an  explicit 
graphical  tool  for  describing  energy  exchange  within  a  system 
and  facilitating  multidisciplinary  exchanges  [11].  To  this  point 
of  view,  bond  graphs  and  VHDF-AMS  share  advantages:  they 
allow  structural  decomposition  of  the  system  into  subsystems 
[12]  and  can  also  use  a  functional  approach;  bond  graph  like 
VHDF-AMS  is  unified  for  all  physical  domains  [13,12]  and  can 
be  easily  connected  to  other  existing  models  such  as  auxiliaries 
like  other  sources  of  energy,  loads  and  power  conversion  devices. 

They  offer  also  a  benefit  on  the  control  design:  they  have 
built-in  causality  assignment,  an  explicit  definition  of  state  vari¬ 
ables  and  a  direct  graphic  interpretation  of  controllability  and 
observability  in  terms  of  causal  paths.  However,  neither  con¬ 
trol  scheme  nor  control  strategies  deduced  from  a  fuel  cell  bond 
graphs  model  have  been  yet  published:  one  can  wonder  if  the 
complexity  of  the  model  itself  does  not  permit  to  find  easily 
control  strategies.  Even  if  bond  graphs  are  a  powerful  tool, 
they  also  have  some  weaknesses.  Firstly,  they  allow  only  energy 
exchanges  which  are  described  by  bonds  which  represent  power 
interaction:  two  variables  effort  and  flux  (equivalent  to  the  vari¬ 
able  across  and  through  in  VHDF-AMS)  are  associated  with 
each  bond.  The  power  is  calculated  by  multiplying  effort  and 
flux:  for  example  in  the  fluidic  domain,  effort  variable  is  the 
pressure  (Pa)  and  flux  is  the  volumetric  flow  rate  (m-3  s).  In 
some  systems,  like  fuel  cells,  it  would  be  more  easy  to  use  mass 
flow  for  mass  conservation  instead  of  volumetric  flow  but  the 
product  between  a  mass  flow  and  a  pressure  does  not  give  a 
power:  therefore,  volumetric  flows  are  used  for  correct  power 
calculation.  Secondly,  bond  graphs,  originally,  do  not  support 
discontinuities:  generally  average  models  of  power  converter  are 
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used  where  switching  is  not  considered  [11]  and  digital  part  of 
the  system  cannot  be  modelled  (e.g.,  FPGA,  digital  controllers). 
Finally,  the  bond  graph  representation,  even  if  it  is  powerful, 
cannot  be  easily  read  by  a  novice  as  contrasted  with  VHDL- 
AMS  where  physical  equations  are  directly  written,  commented 
and  possible  to  be  understood  by  a  novice. 

2.2. 7.2.  Electrical  equivalent  circuits.  In  an  attempt  to  support 
the  modelling  and  simulation  of  nonelectrical  as  well,  several 
modelling  methods  using  energy  equivalences  between  the  elec¬ 
trical  domain  and  other  domains,  such  as  mechanical,  thermal, 
or  fluidic  domains,  have  been  adopted.  The  basic  method  for 
creating  a  new  model  in  a  Spice-like  environment  is  to  combine 
pre-existing  building  blocks  to  achieve  the  desired  functionality: 
voltage  source,  current  controlled  voltage  sources  [14],  resistor, 
capacitance  and  even  semi-conductors  [15]. 

The  electrical  equivalent  fuel  cell  models  can  be  divided  in 
two  parts:  the  fuel  cell  electrical  part  model  which  contains 
only  equivalent  voltage  sources,  current-controlled  voltages  and 
resistor  for  activation,  ohmic  and  concentration  losses.  The  dou¬ 
ble  layer  charging  effect  is  modelled  by  a  capacitance.  The  most 
well-known  simple  equivalent  circuit  model  have  been  presented 
and  explained  in  many  papers  and  books  [16-19].  Other  elec¬ 
trical  equivalent  circuit  models  have  been  presented  and  are 
summarised  in  [20] :  most  of  them  are  mainly  based  on  empiri¬ 
cal  results  where  parameters  are  identified  by  current  interrupt 
technique  or  impedance  spectrometry.  However,  empirical  rela¬ 
tionships  do  not  provide  an  adequate  physical  understanding  of 
the  phenomena  inside  the  cell. 

Fuel  cell  gas  channels  have  been  modelled  using  equivalent 
electrical  circuits.  Mole  or  mass  conservation  of  each  species 
(H2,  H2O,  O2)  are  modelled  using  three  different  electrical  cir¬ 
cuit  where  the  voltage  sources  represent  pressures,  the  current 
sources  represent  the  flows  [14,21].  Other  components  such 
as  resistors  and  capacitor  represents  the  fluidic  resistance  and 
capacitance  of  the  channels.  The  model  of  Famouri  and  Gem- 
men  [14]  includes  flow  transients,  reactant  partial  pressures  and 
loss  mechanisms  within  the  fuel  cell.  Hernandez  and  Diong  [21] 
takes  into  account  some  nonlinearities  like  the  vapour  saturation 
pressure  but  results  show  a  divergent  behaviour  if  the  operat¬ 
ing  point  zone  migrates  far  away  from  the  conditions  used  for 
model  parameters  identification  because  capacitances  and  resis¬ 
tances  which  are  used  in  the  model  are  constant.  Wang  et  al.  [22] 
developed  also  a  complex  dynamic  fuel  cell  model  using  elec¬ 
trical  analogies  taking  into  account  double  layer  charging  effect 
and  the  thermodynamic  characteristic  inside  the  fuel  cell.  The 
Wang  model  has  been  successfully  used  in  a  distributed  gen¬ 
eration  system  model  including  power  converters  and  control 
[23]. 

The  electrical  equivalent  circuits  have  a  greater  limitation 
because  they  needs  pre-defined  buildings  blocks.  But  what  if  the 
designer  needs  special  effects  for  these  models  like  temperature 
dependency?  How  about  hydraulic  or  mechanical  components 
that  cannot  be  easily  described  with  standard  building  block 
components?  Designers  may  have  several  requirements  that  are 
possible  to  model  using  pre-defined  building  blocks.  Moreover, 
while  the  function  can  be  realised  in  this  manner,  it  is  rather  hard 


to  parameterise  and  reconfigure  if  needed:  working  with  this  type 
of  macromodel  can  be  a  real  limitation  to  system  designer  who 
needs  a  model  quickly. 

2.2. 1.3.  Signal  flows.  Signal  flows  modelling  is  the  most  used 
approach  in  control-oriented  models:  many  fuel  cell  models  have 
been  published  using  Matlab/Simulink  environment.  Lukas  et  al. 
[24]  developed  a  control-oriented  dynamic  simulation  model  of  a 
direct  reforming  molten  carbonate  fuel  cell  (MCFC)  powerplant 
and  extended  further  the  model  with  experimental  validation  for 
several  load  points  in  [25].  Uzunoglu  and  Alam  [26]  developed 
a  hybrid  fuel  cell  and  supercapacitors  dynamic  model  where  the 
fuel  cell  model  takes  into  account  a  methanol  reformer  (modelled 
by  a  transfer  function)  and  an  electrical  fuel  cell  steady  state 
characteristic;  this  approach  is  a  functional  approach  which  is 
well  suitable  to  implements  controllers  (proportional-integral 
controller  in  this  case)  to  control  the  overall  system.  However, 
in  this  kind  of  models,  there  is,  on  the  one  hand,  the  graphical 
part  (Simulink)  of  the  model  (see  [26]  [Fig.  2])  and,  on  the  other 
hand,  a  file  which  contains  all  the  numerical  parameters.  This 
makes  the  model  really  hard  to  read  and  it  cannot  be  used  easily 
by  someone  who  wants  to  simply  re-parametrise  or  improve 
some  parts  of  the  model. 

Pasricha  and  Shaw  [27]  extended  the  static  current  voltage 
description  to  include  temperature  dependence  with  a  dynamic 
modelling  of  the  membrane  temperature.  Correa  et  al.  [28]  devel¬ 
oped  a  computer-controlled  high-power  converter  to  emulate 
static  and  dynamic  fuel  cell  characteristics:  the  model  includes 
membrane  temperature  and  humidity  dependencies,  efficiency, 
reactant  flows,  cooling  air  fans  and  water  pump.  This  model 
has  been  implemented  in  a  Lab  View  (National  instruments) 
environment  which  allows  signal-flow  graphical  programming. 

Stephanopoulou  and  co-workers  [29-32]  developed  one  of 
the  most  advanced  control-oriented  nonlinear  fuel  cell  model 
until  now.  The  model  takes  into  account  dynamic  fuel  cell 
reactant  supply  based  on  lumped- volume  filling  dynamics,  aux¬ 
iliaries’  components  (compressor,  manifold,  static  air  cooler, 
static  humidifier).  The  fuel  cell  stack  model  is  based  on  a  func¬ 
tional  approach  and  it  is  separated  into  four  functionalities:  the 
stack  voltage,  the  cathode  flow,  the  anode  flow  and  the  mem¬ 
brane  hydratation.  This  kind  of  approach  is  very  powerful  and 
it  is  well  suitable  for  control  design  and  analysis.  But  it  also  has 
some  drawbacks  because  the  different  functionalities  are  highly 
inter-dependents  and  they  cannot  be  easily  simulated  separately 
because  of  feedbacks  loops.  One  can  notice  that  in  the  functional 
approach  each  functionality  is  a  physical  domain:  humidity,  elec¬ 
trical,  fluidic  functionalities  are  separated.  Moreover,  if  a  system 
designer  wants  to  improve  one  part  of  the  model,  for  example  the 
cathode  model,  to  take  into  account  a  new  phenomena,  he  has 
to  modify  not  only  the  cathode  flow  functionality  (gas  channels, 
supply  manifold,  return  manifold),  but  also  the  voltage  function¬ 
ality.  In  this  case,  how  can  he  knows  which  part  of  the  model 
he  has  to  modify  or  not?  Which  part  of  the  voltage  function 
corresponds  to  the  cathode  overpotential?  How  can  he  be  sure 
he  modified  all  what  had  to  be  modified?  How  can  he  test  and 
prove,  independently  of  the  other  parts,  the  validity  of  the  new 
phenomena  he  took  into  account? 
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Fig.  6.  Fuel  cell  structure  abstraction  levels:  (a)  stack  level;  (b)  layers  level. 


In  the  structural  approach,  as  it  will  be  presented  in  the  next 
part,  these  kind  of  problems  do  not  arise:  if  the  system  designer 
wants  to  change  the  cathode  model,  he  knows  that  all  the  multi- 
domain  phenomena  (electrical  voltage,  fluidic,  thermal,  etc.)  of 
the  cathode  are  taken  into  account  in  the  cathode  component. 

2.2.2.  Top-down  approach 

The  modelling  approach  can  be  described  by  Fig.  6.  Two  level 
of  abstraction  are  given:  the  first  one  Fig.  6(a)  can  be  called  the 
’’stack  abstraction  level”  where  the  fuel  cell  stack  is  seen  only  as 
a  voltage  source:  in  this  case  all  the  fuel  cell  voltage  is  computed 
in  the  same  function  taking  into  account  the  losses  as  it  is  done 
in  many  papers. 

The  second  level  (Fig.  6(b))  is  more  refined  and  can  be  called 
“layer  abstraction  level”  because  it  takes  into  account  more  parts 
and  refined  phenomena.  In  this  article,  this  ’’layer  abstraction 
level”  is  used.  Several  structural  parts  are  distinguished: 

•  The  cathode  gas  channels:  it  is  a  simple  pipe  where  gas  can 
enter  and  leave.  In  this  model,  pressure  losses  are  not  taken 
into  account  and  can  be  included  in  a  more  refined  model; 

•  The  cathode  gas  diffusion  layer  (GDL)  where  water  and  oxy¬ 
gen  diffuse; 

•  The  cathode  throttle  (or  valve)  which  permits  to  regulate  the 
cathode  pressure; 


•  The  cathode  catalyst  layer  where  the  oxygen  is  consumed  and 
the  water  is  produced; 

•  The  membrane  where  water  transport  phenomena  are  taken 
into  account; 

•  The  anode  catalyst  layer  where  hydrogen  is  consumed; 

•  The  anode  gas  diffusion  layer  where  hydrogen  and  water 
diffuse; 

•  The  anode  gas  channels. 

Each  component  imposes  a  voltage  which  can  be  dependent 
on  any  quantity  of  the  component  (e.g.,  current,  partial  pressure) 
and  the  overall  stack  voltage  is  obtained  by  plugging  the  com¬ 
ponents  in  series.  Each  component  can  be  modelled,  described, 
tested  and  validated  independently . 

2.2.2. 1.  Assumptions.  In  the  first  modelling  steps,  some 
assumptions  are  given.  Building  a  more  refined  model  will  con¬ 
sist  to  remove  one  by  one  the  following  assumptions: 

HI.  Diffusion  is  considered  in  steady  state. 

H2.  Linear  diffusion  (concentration  gradient). 

H3.  No  nitrogen  diffusion. 

H4.  No  pressure  losses  in  the  cathode  gas  channels. 

H5.  Hydrogen  pressure  is  constant  and  perfectly  regulated. 


440 


B.  Blunier,  A.  Miraoui  /  Journal  of  Power  Sources  177  (2008)  434^150 


H6.  Temperature  is  homogeneous  in  the  entire  stack. 

H7.  Cathode  volume  is  constant. 

H8.  No  water  accumulation  in  the  anode  and  cathode  catalyst 
layers. 

H9.  Anode  catalytic  activity  is  considered  so  high  that  anode 
overpotential  is  neglected.  In  this  case,  all  activation  losses  occur 
at  the  cathode  catalyst  layer. 

H10.  Even  more  or  less  important  disparities  exist  between 
the  cells  of  a  stack,  it  is  assumed  here  that  n- cell  stack  is  n 
times  the  behaviour  of  the  equivalent  mean  cell.  Note:  there 
are  functionalities  in  VHDL-AMS  (generate  function)  which 
could  permit  to  instantiate  n  cells  to  build  an  overall  stack  but 
this  is  not  investigated  here. 

Hll.  The  gases  are  assumed  to  be  perfect. 

H12.  Water  does  not  leave  the  stack  in  the  liquid  form  but  only 
in  vapour  form:  liquid  water  accumulates  or  evaporate  in  the 
cathode  gas  channels  [30] . 

According  to  H8  and  HI 2,  the  model  cannot  capture  electrode 
flooding  phenomenon.  The  effect  of  excess  of  liquid  water  on 
oxygen  mass  transport  rates  within  the  cathode  is  really  impor¬ 
tant.  Two-phase  flow  model  for  PEMFCs  is  currently  an  area  of 
active  research.  By  neglecting  this  phenomenon  we  introduce 
significant  errors  in  the  results  of  the  PEMFC  cathode  water 
distribution.  It  also  affects  the  cathode  overvoltage  computa¬ 
tion,  making  the  model  less  realistic  in  given  conditions.  Most 
of  the  two-phase  models  are  solved  using  CFD  tools  which  solve 
spatial  partial  differential  equations  (e.g.,  Naviers-Stokes).  An 
analytical  model  has  been  given  by  Baschuk  [33]  and  can  be 
added  in  a  more  refined  model. 

However,  the  tests  have  been  performed  on  a  low  power  fuel 
cell  (1.2  kW)  which  works  at  nearly  atmospheric  pressure  and 
a  high  stoichiometry  ratio  (the  measured  ratio  is  at  least  5).  In 
these  conditions  there  are  very  few  chances  to  flood  electrodes 
of  the  fuel  cell:  the  model  has  been  validated  in  these  conditions. 
In  other  conditions,  for  example  at  high  current  and  low  speed 
air  stream,  the  model  will  not  be  able  to  predict  water  flooding. 


2.2.3.  Gas  channels 

2.2.3. 1.  Cathode.  According  to  hypothesis  H4,  the  model  does 
not  take  into  account  pressure  losses  in  the  cathode  channels. 
The  total  pressures  at  the  inlet  and  outlet  manifold  (superscript 
m)  pfnm  and  respectively,  are  equal  to  the  total  cathode 
pressure  pc\ 

#4  =►  P?nm  =  =  Pc  =  Pc0l + pg2 + /40  (i) 

where  Pq ,  p^2  and  P§2q  are  the  cathode  oxygen,  nitrogen  and 
vapour  manifold  partial  pressures,  respectively. 

Considering  all  the  gases  as  perfect  gases  (Hll),  it  can  be 
written, 


niRT 


where  pi  and  ni  are  the  partial  pressure  and  number  of  moles  of 
the  species  i,  respectively.  cathode  channels  volume  and  T  is 
the  temperature. 

Derivating  Eq.  (2)  gives, 

d pi  RT  dnt  Rnt&T  niRT  dl/^ 

1  — _ 1  4 - i _ _ _  ni 

d  t  V0  d  t  V0  d  t  yC2  d  t 

According  to  H7  (dV^/dt  =  0),  the  previous  equation  can  be 
simplified: 

d  pi  _  RT^  d Tf  Rn^  d T 
d t  V0  dt  d t 

d pi  RT  1  d mi  Rmi  dT 

1  = - -  -| - - -  (5) 

dt  Mi  dt  Mi  dt 

where  mz-  and  Mz  are  the  mass  and  molar  mass  of  the  species  i, 
respectively. 

Each  species  inside  the  cathode  channels  can  enter  inside  the 
channels  (subscript  in),  go  outside  the  fuel  cell  (subscript  out)  or 
go  inside  the  catalyst  layer  (subscript  inside).  The  mass  balance 
for  each  species  can  be  written, 

dm/ 

=  Qi, in  H”  Qi, out  H”  Qi, inside  with  i  £  {O2,  N2,  H2O}  (6) 
dt 

It  has  to  be  noticed  that  there  are  no  assumption  on  the  sign 
of  the  mass  flows:  the  mass  flows  are  imposed  by  (1)  external 
devices  like  the  compressor  and  the  humidifier  (2)  the  valve, 
(3)  the  catalyst  layer  and  they  can  be  either  positive  or  negative 
depending  on  the  operating  point. 

Eq.  (5)  applies  for  gases  only  whereas  (6)  applies  for  gases 
and  liquids.  The  case  where  water  is  in  liquid  form  has  to  be 
distinguished  in  order  to  take  into  account  saturation  of  vapour 
in  the  air  (Fig.  7). 
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Fig.  7.  Cathode  gas  channels. 
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Two  equations  have  to  be  added  [30].  The  first  one  (7)  checks 
if  the  total  mass  of  water  inside  the  control  volume  (gas  channels) 
is  greater  than  the  mass  of  vapour  which  can  saturate  the  air 
m^’2Q  .  If  the  total  mass  of  vapour  is  lower  than  m^Q  ,  then 
all  the  water  is  in  the  vapour  form;  in  the  other  case,  the  air  is 
saturated  (jbh2q  =  ^n2o  )■ 

The  second  equation  (8)  links  the  total  mass  of  water,  the 
vapour  (subscript  V)  and  liquid  (subscript  C)  mass  water. 

C,V  _  fmH20  mH20  -  mH2Oat  (1, 

H2°  1  mH2oat  elsewhere 

mH2o  =  mH20  +  mH2o  (^) 


The  mass  of  water  which  saturates  the  air  m^{2Q  is  calculated 
as  follows, 


C,  Sat  _  )  vc  Mh2o 

mH20  r  j> 


(9) 


From  the  oxygen  partial  pressure,  the  fuel  cell  EMF  cathode 
contribution  can  be  computed  [35]: 


Ec  =  n 


1.229  -0.85  x  10“3(T 


298.15) 


RT 

+ - In 

IF 


(17) 


2.23.2.  Anode.  The  model  of  the  anode  gas  channels  is  based 
on  the  cathode  gas  channels  model  but  according  to  H5,  it  is 
assumed  that  the  hydrogen  partial  pressure  is  constant.  The  fuel 
cell  is  in  dead-end  mode.  Nor  hydrogen,  nor  water  can  leave  the 
fuel  cell  at  the  anode  gas  channels. 

Considering  this,  it  can  be  written: 


A,  V  _  f  mH20  if  mH20  —  mU20 

H2°  1  elsewhere 


(18) 


where  /?sat  is  the  vapour  saturation  pressure  given  by  [34], 


A  A, V  ,  A ,C 

mH20  —  mH20  +  mH20 


(19) 


,  (  ps*t(T)  \ 

0§1°  (l.O  X  io5 ) 

=  -2.1794  +  0.02953  (T  -  273.15) 

-  9.1837  x  1(T5  (T  -  273. 15)2 

+  1.4454  x  10“7  (T  -  273. 15)3  (10) 

where  temperature  T  and  pressure  psat  are  in  Kelvin  and  Pascal, 
respectively. 

Finally  a  set  of  differential  equations  is  obtained: 

dpg2  _RT  1  dmg2  |  Rm^  dT 
d  t  Mq2  d  t  \f  Mq2  d  t  1 


The  mass  of  water  which  saturates  the  air  m^Qlt 
as  follows, 

A,  Sat  _  PSaxOTAMiW 

mH20  r  j 

dPH2o  _  RT  1  dmH2g  Rtn^f  d T 
d  t  Mh2o  d  t  Mh2o  d  t 

d  jAm 

2  =  o 

d  t 

dmH2Q  _  A 

1  f  ^H20, inside 


is  calculated 


(20) 

(21) 

(22) 

(23) 


d p^2  _  RT  1  dmg2  Rm^2  d T 
d  t  Mn2  d  t  Mn2  d  t 

d  I}<h2o  _  R 1  \  d  Wh2o  U oth2o  d  T 

d  t  )/-  Mh2o  d  t  V0  Mh2o  d  t 

dwg2  =ac  ,  c  ,  c 

W2,in  '  ^02,out  '  inside 

dmN2  _  C  C 

7  —  *7  N  2 .  i  n  +  7N2,0Ut 


(12) 


(13) 

(14) 

(15) 


dwH2Q 

d  t 


c  c  c 

^H20,in  "f"  te20,out  "f"  fe20, inside 


(16) 


According  to  H3,  no  nitrogen  diffuses  from  the  channels  to  the 
catalyst  layer  that  is  why,  in  the  mass  conservation  equation  of 
nitrogen  (15),  inside  is  not  taken  into  account. 


where  VA  is  the  anode  volume.  According  to  H5  the  hydrogen 
inlet  mass  flow  #H2?n  cann°t  be  controlled:  it  is  imposed  by  the 
anode  catalyst  layer  downstream  of  the  anode  gas  channels;  on 
the  contrary,  the  inlet  air  mass  flow  in  the  cathode  gas  channels 
is  imposed  by  the  compressor  and  humidifier  upstream  of  the 
gas  channels. 

From  the  hydrogen  partial  pressure,  the  fuel  cell  EMF  anode 
contribution  can  be  computed  [35]: 


a  RT 

Ea  =  n -  In 

2  F 


Pu2 


101,325 


(24) 


The  component  of  the  anode  gas  channels  is  given  in  Fig.  8. 


2.2.4.  Control  valve  or  outlet  throttle 

Depending  on  the  operating  mode,  the  pressure  inside  the 
cathode  can  be  controlled  or  not:  the  throttle  opening  area  can 
be  set  to  a  constant  or  can  be  used  as  an  extra  control  variable 
to  regulate  the  cathode  pressure  [30] . 


442 


B.  Blunier,  A.  Miraoui  /  Journal  of  Power  Sources  177  (2008)  434^150 


ea 


*7h2  -inside 


Ph2 


^H2  0 ’inside 


Ph2o 


9H2-in 


?H20.iri 


Ph2o 


Ph2 


Fig.  8.  Anode  gas  channels. 


KN2  = 


_ Xn2  MNi _ 

X02  Mo2  +  xn2  Mn2  +  xh20  ^h2o 


Kh2o  = 


_ Xh2o  ^h2o _ 

X02  Mo2  +  Xn2  Mn2  +  XH20  Mh2o 


(28b) 

(28c) 


where  xo2>  Xn2  and  xh2o  are  the  mole  fraction  of  oxygen, 
nitrogen  and  vapour,  respectively. 

The  molar  fractions  computation  changes  if  the  mass  flow  is 
positive  in  >  Ptot.om)  or  negative  (p^  in  <  out)evenif 
in  a  fuel  cell  the  mass  flow  should  always  be  positive  (from  the 
cathode  to  the  atmosphere),  this  case  is  taken  into  account: 


Xo2  = 


Xn2  = 


PO2,  in/ ^tot,  in 

if 

^tohin  —  Ptot,  out 

V  /V 

P 02,  out/ Ptot,  out 

elsewhere 

PN2,  in/ Ptot,  in 

if 

Ptot,  in  —  Ptot,  out 

V  /V 

^N2,  out/ P tot,  out 

elsewhere 

(29a) 

(29b) 


The  valve  or  the  throttle  (see  Fig.  9)  can  be  modelled  by  a 
linearised  form  of  the  subcritical  nozzle  flow  equation  as  given 
in  [30],  that  is, 

^tot  —  ^  (Ptot,  in  —  Ptot,  out)  (25) 


XH20  = 


Pu20,in/Plt,m  if  Ptot,  in  >  Ptot,  out 
Ph2o,  out/ Ptot,  out  elsewhere 


(29c) 


The  total  pressure  at  the  inlet  and  outlet  depend  also  on  the  three 
species: 


where  g^t  is  the  total  mass  flow  rate  going  through  the  nozzle 
and  pVQt  in  and  p^ot  out  the  inlet  and  outlet  total  pressures  across 
the  nozzle. 

Three  species  (O2,  N2  and  H2O)  are  inside  the  gas  and  the 
total  mass  flow  can  be  written, 


V  _  V  .  V  .  V 
tftot  —  TO2  '  #N2  '  ^H20 

(26) 

With, 

V  V 

^o2  —  yo2  Qtot 

(27a) 

V  V 

<?N2  —  KN2  tftot 

(27b) 

#H2Q  =  KH20  <liot 

(27c) 

where  yo2,  y^2  and  xh2o  are  the  mass  fraction  of  oxygen,  nitro¬ 
gen  and  vapour,  respectively. 

The  mass  fraction  of  the  three  species  are  computed  as  fol¬ 
lows: 

xo2  mo2 

yo2  — 

Xo2  Mo  >  +  Xn2  mn2  +  Xh2o  Mh2o 

(28a) 

v  v  .  v  .  v  /ork  x 

Plot,  in  -  Po2,  in  +  ^N2,  in  +  ^H20,  in  (30a) 

Ptot,  out  PO2,  out  P^2,  out  PH2O,  out  (30b) 


2.2.5.  Catalyst  layers 

Catalyst  layers  impose  the  mass  flow  rates  in  function  of  the 
stack  current  (/).  It  is  assumed  that  no  pressure  losses  occur 
in  the  catalyst  layers:  the  inlet  and  outlet  vapour  pressures  are 
equal.  The  component  of  a  catalyst  layer  is  given  in  Fig.  10. 


2. 2. 5.1.  Cathode  catalyst  layer.  The  electrical  current  imposes 

an  oxygen  mass  flow  to  the  cathode  gas  diffusion  layer  (subscript 

C  c 

GDL).  The  oxygen  mass  flow  ^q2gdl  taken  from  the  cathode 
GDL  to  the  cathode  catalyst  layer  is  given  by 


C,  c 

4o2,gdl 


Mq2 


In 
4 F 


(31) 


Fig.  9.  Valve  or  throttle. 


Fig.  10.  Catalyst  layer. 
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At  the  same  time,  water  is  produced  inside  the  catalyst  layer, 

A O,  prod  =  MH20  ^  (32) 

Water  vapour  mass  flow  coming  from  the  membrane  to  the  cat- 
alyst  layer  gH’2o  m  *s  not  known  a  priori ,  it  is  imposed  by  the 
membrane  and  will  be  introduced  in  the  membrane  model  (see 
Section  2.2.9). 

Moreover,  according  to  H8  no  water  accumulates  in  the  cat¬ 
alyst  layer  (dmj^o/d  t  =  0). 

Water  mass  balance  can  therefore  be  written, 


2.2.6.  Gas  diffusion  layers 

The  relation  between  the  current  density  j  and  the  reactants 
molar  flux  J  (mol  cm-2  s-1)  coming  from  the  channels  to  the 
catalyst  is 


(39) 


where  ne  is  the  number  of  electrons  entering  in  the  reaction. 

In  steady  state  (HI)  the  reactant  molar  flux  (from  channel  to 
the  catalyst)  is  equal  to  the  product  molar  flux  (from  the  catalyst 
layer  to  the  channels). 


C,  c  _  /  C,  c  .  C,  c  \  /oo\ 

^H20,GDL  —  -b?H20,M  +  <^0,  prod! 

The  minus  sign  is  to  take  into  account  mass  flow  conventions.  If 

C  c 

the  mass  flow  of  water  produced  gH2o  prod  at  catalyst  layer 
is  always  positive,  the  water  mass  flow  exchanged  between  the 
membrane  and  the  catalyst  q^o  M  can  edher  positive  (water 
comes  from  the  membrane)  or  negative  (water  is  absorbed  by 
the  membrane). 

Activation  losses  occur  in  the  catalyst  layer  and  are  computed 
as  follows: 


vif  =  -nA  In  I  J- 


(34) 


where  A  and  b  are  empirical  coefficients.  The  current  density  j 
is  defined  as  follows, 


J  = 


n  I 
3tot 


(35) 


where  Stot  is  the  total  stack  active  area. 


2.2.6. 1.  Diffusion.  The  Fick  law  describes  the  diffusion  as  fol¬ 
lows  [19]: 


dc 

J  =  ~D—  (40) 

ax 

where  J  is  the  diffusion  molar  flux,  c  the  molar  concentration 
and  D  is  the  diffusion  coefficient. 

According  to  H2,  Eq.  (40)  can  be  written, 


/=  Dci(S)-Ci(  Q) 


(41) 


where  q(x )  is  the  reactant  concentration  of  i  at  abscissa  x  and  8 
is  the  gas  diffusion  layer  thickness. 

It  has  to  be  noticed  that  J  can  be  either  positive  (reactant 
species)  or  negative  (product  species)  as  shown  in  Fig.  11. 

The  flux  of  the  species  /,  //,  can  be  computed  from  the  mass 
flow  qi  according  to: 


Ji  —  qi  Mi  Stot 


(42) 


2. 2. 5. 2.  Anode  catalyst  layer.  The  electrical  current  imposes  a 

hydrogen  mass  flow  to  the  anode  gas  diffusion  layer  (subscript 

A  c 

GDL).  The  hydrogen  mass  flow  #h^gdl  taken  from  the  anode 
GDL  to  the  anode  catalyst  layer  is  given  by, 


<?h^gdl  —  ^h2 


In 

2F 


(36) 


No  water  is  produced  in  the  anode  catalyst  layer.  Water  vapour 
mass  flow  coming  from  the  membrane  to  the  catalyst  layer 
<?h’2o  m  is  not  known  a  priori ,  it  is  imposed  by  the  membrane 
itself  and  will  be  introduced  in  the  membrane  model  (see  Section 
2.2.9). 

Moreover,  according  to  H8  no  water  accumulates  in  the  cat¬ 
alyst  layer  (dmH2o/dt  =  0). 

Water  mass  balance  can  therefore  be  written, 


A.  c  C  c 

^H20,  GDL  =  -<?H20,M 


(37) 


The  minus  sign  is  to  take  into  account  mass  flow  conventions: 
the  water  mass  flow  exchanged  between  the  membrane  and  the 
anode  catalyst  qff2o  M  can  be  either  positive  (water  comes  from 
the  membrane)  or  negative  (water  is  absorbed  by  the  membrane). 

According  H9,  there  are  no  voltage  losses  in  the  anode  cata¬ 
lyst  layer: 


vtf  =  o 


(38) 


K»ncDL  with  i  €  {CtA} 


reactant 


fv(o) 

r 

Gas  channel 
interface  product 


Pp(  0) 


Cr(0)  9r>0 


Gas  Diffusion 
Layer  (GDL) 

Cp(0)  %  <  0  Cp(8) 


Pr(S) 


Catalyst 

interface 


pP(<s) 


fusion  layer. 
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Table  1 

Critical  properties  of  gases 


Species 

Molar  mass 

Tc  (K) 

Pc  (atm) 

h2 

2.016 

33.3 

12.80 

Air 

28.964 

132.4 

37.0 

n2 

28.013 

126.2 

33.5 

02 

31.999 

154.4 

49.7 

h2o 

18.015 

647.3 

217.5 

where  Mz-  is  the  molar  mass  of  the  species  i  and  Stot  is  the  total 
active  area  of  the  stack. 

Concentration  profiles  in  steady  state  in  the  gas  diffusion 
layer  are  given  in  Fig.  1 1 

Concentration  of  the  species  i  along  the  abscissa  x  can  also 
be  written  as  a  function  of  the  molar  fraction  Xi(x)  and  the  total 
pressure  p  or  partial  pressure  Pi(x), 

Ci{x)  =  -Jj;  Xi(x)  =  d-  Pi(x)  (43) 

For  low  pressures,  the  binary  diffusion  coefficient  D;y  (two 
species  i  and  j)  depends  on  the  pressure  and  the  temperature 
according  to  the  following  Slattery  and  Bird  gas  law  [36]: 


pDij  =  a 


( PciPcj)l/3{TciTcj)5/l 2 


X 


_L  j_y/2 

Mi  +  Mj) 


(44) 


where  p  is  the  total  pressure  (atm),  the  binary  coefficient 
(cm2 */s)  and  T  is  the  temperature  (K).  Mx  is  the  molar  mass  of  the 
species  x,  Tcx  and  pcx  are  the  critical  temperature  and  pressure, 
respectively,  of  the  species  x  (x  e  {/,  j }).  The  coefficients  a  and 
b  differ  if  one  of  the  two  species  is  a  polar  gas  or  not  (water  in 
this  case).  Then, 

a  =  2.745  x  1(T4  and  b  =  1.823  (45) 


if  the  pair  of  gas  contains  nonpolar  gases  and 
a  =  3.640  x  1(T4  and  b  =  2.334  (46) 


if  the  pair  of  gas  contains  a  polar  gas  (H2O). 

The  critical  temperatures  and  pressures  of  the  different 
species  interacting  in  a  fuel  cell  are  given  in  Table  1 . 

For  porous  media,  the  binary  diffusion  coefficient  has  to  be 
corrected  to  account  for  the  effects  of  the  pore  walls  (poros¬ 
ity,  tortuosity).  Usually,  this  is  accomplished  by  employing  a 
modified  or  effective  diffusivity,  also  known  as  the  Bruggemann 
correction  for  the  diffusion  coefficients  [34,19]: 


K  = 


or  Dij-  for  high  temperatures 


(47) 


fuel  cell  electrodes)  which  describes  the  additional  impedance 
do  diffusion  caused  by  tortuous  path  flow. 


2.2.7.  Anode  gas  diffusion  layer  model 

Two  species,  H2  and  H2O,  are  present  at  the  anode  gas  dif¬ 
fusion  layer  interface  (hypothesis  H3  permits  to  neglect  the 
nitrogen  diffusion)  and  Eq.  (47)  is  used  to  compute  H2q. 

According  to  Eqs.  (41)  and  (43),  the  partial  pressures  of 
the  hydrogen  and  the  vapour  at  the  catalyst  interface  can  be 
expressed: 

,  jA  RT  , 

Ph2(Sa)  =  ph2(0)  -  (48) 

^H2,H20 

jA 

Ph2o($A)  =  /?h2o(0)  -  - 8a  (49) 

^h2,h2o 

where  the  binary  diffusion  coefficient  D^2  H2q  is  calculated 
from  Eqs.  (44),  (45)  and  (47). 

The  molar  flux  J^2  and  J^2q  are  computed  from  the  cor¬ 
responding  mass  flows  q^2  and  respectively,  with  Eq. 

(42). 

The  anode  concentration  voltage  drop  is  deduced  from  the 
catalyst  and  bulk  concentrations  (or  partial  pressures): 


VA,  GDL 
*conc 


— n  C  In 


(  Ph2(0)  \ 
Uh2(<$a)) 


where  C  is  an  empirical  coefficient. 


(50) 


2.2.8.  Cathode  gas  diffusion  layer  model 

Two  species,  O2  and  H2O  are  present  at  the  cathode  gas 
diffusion  layer  interface  (hypothesis  H3  permits  to  neglect  the 
nitrogen  diffusion)  and  Eq.  (47)  is  used  to  compute  Dq2  H2q. 

According  to  Eqs.  (41)  and  (43),  the  partial  pressures  of 
oxygen  and  vapour  at  the  catalyst  interface  can  be  expressed: 

j  c  y 

Po2(Sc)  =  po2(0)  -  -% - <5C  (51) 

^02,H20 

j c  y  p 

pn2o(f) = 7>h2o(0)  -  yy —  ^  (52) 

do2  h2o 

where  the  binary  diffusion  coefficient  Dq2  H2q  is  calculated 
from  Eqs.  (44),  (45)  and  (47). 

The  cathode  concentration  voltage  drop  is  deduced  from  the 
catalyst  and  bulk  concentrations  (or  partial  pressures): 

t/C, GDL  P° 2(0)  (C2) 

<53) 


where  6  stand  for  the  porosity  of  the  structure  (around  0.4  for 
fuel  cell  electrodes1)  and  r  is  the  tortuosity  (from  1.5  to  10  for 


1  A  porosity  of  0.4  means  that  40%  of  the  total  electrode  volume  is  occupied 

by  pores.  In  open  space,  the  porosity  is  1  and  the  normal  diffusion  coefficient  is 

employed. 


2.2.9.  Membrane 

The  membrane  model  (see  Fig.  12)  developed  here  is  the 
Springer  isothermal  and  one-dimensional  model  given  in  [34] 
where  the  authors  empirically  determine  relationships  correlat¬ 
ing  membrane  conductivity  and  electrode  porosity  with  water 
content  in  the  Nafion  membrane. 
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The  membrane  resistance  is  deduced  from  its  conductivity. 
Since  membrane  conductivity  is  highly  dependent  on  its  water 
content,  it  is  essential  to  know  how  water  content  varies  across  a 
Nation  membrane.  Two  opposite  phenomena  occur  in  the  mem¬ 
brane  and  have  to  be  taken  into  account:  the  electro-osmotic  drag 
and  the  back  diffusion  of  water. 


2.2.9. 7.  Water  content  of  the  membrane.  The  water  content  A 
in  a  Nation  membrane  is  not  uniform  as  it  will  be  shown  in  next 
parts.  It  is  defined  as  the  ratio  of  the  number  of  water  molecules 
to  the  number  of  charged  sites.  Experimental  results  show  that  A 
can  varies  from  0  (completely  dehydrated  membrane)  to  22  (full 
saturation,  under  certain  conditions)  [19].  In  order  to  compute 
the  water  profile  in  the  membrane,  boundary  conditions  have  to 
be  known.  The  water  content  is  related  to  the  water  activity  on 
the  faces  of  the  membranes  according  to  the  following  law  [34] : 


f  0.043  +  17.81  <2h2o  —  39.85  a^20  36.0  ^^j2q  if  0  c  ^h2o  —  1 

\  14  +  1.4(<2h2o  —  1)  if  1  <  au2o  <  3 

(54) 


is  quantified  by  the  electro-osmotic  drag  coefficient  ft  drag*  which 
is  defined  as  the  number  of  water  molecules  accompanying  the 
movement  of  each  proton  [19].  It  is  commonly  assumed  that 
n drag  varies  linearly  with  A  as  [34]: 

«drag  =  «drag  f  tor  0  <  A  <  22  (57) 

where  ~  2.5  is  the  electro-osmotic  drag  coefficient  in  fully 
hydrated  Nation. 

Other  drag  coefficients  expressions  have  been  reported  in  the 
literature  [37,38]  but  Eq.  (57)  have  been  widely  used  [30,34] 
and  have  shown  good  results. 

The  water  drag  flux  /^0  drag  (mol  cm-2  s-1)  from  anode 
to  cathode  as  a  function  of  the  current  density  j  (A  cm-2)  and 
electro-osmotic  drag  coefficient  ft  drag  *s  given  by 

^0,drag  =  2^drag  (58) 


2.2.93.  Back  diffusion  of  water.  When  the  concentration  of 
water  at  the  cathode  side  is  higher  than  the  one  at  the  anode 
side,  water  diffuses  back  from  the  cathode  to  the  anode.  This 
phenomenon  counterbalances  the  effect  of  the  electro-osmotic 
drag.  The  water  back-diffusion  flux  can  be  determined  by  [34]: 


tM 

JH20,  diff 


Pdry  ~  d  A 

D  A 

Mm  ax 


(59) 


where  pdry  is  the  Nation  dry  density  (pdry  = 
0.00197  kg/m3  [19]),  Mm  is  the  Nation  equivalent  weight 
(Mm  =  1.0  kg/mol  [19])  and  v  is  the  direction  through  the 
membrane  thickness  from  the  anode  to  the  cathode. 

The  water  diffusivity  D\  (cm2/s)  is  not  constant;  it  depends 
on  the  temperature  T  and  water  content  A  of  the  membrane 
according  the  following  law  [34] : 


These  values  are  given  for  a  temperature  of  30  °  C  but  it  has  been 
assumed  in  [34]  that  it  is  valid  also  up  to  80  °C. 

The  water  activity  is  calculated  from  the  partial  water  pressure 
Ph2o  and  saturation  pressure  psat- 


ftH20 


Ph2o 

P$at(T) 


(55) 


where  the  vapour  saturation  pressure  psat  is  given  by  Eq.  (10). 

From  the  boundary  water  partial  pressure  conditions  at  the 
cathode  (superscript  C)  and  anode  (superscript  A),  the  average 
Nation  water  content  can  be  found: 


(A)  = 


AC  +  A^ 


(56) 


where  Ac  is  found  from  the  water  partial  pressure  at  the  cathode 
GDL  side  and  Aa  from  the  water  partial  pressure  at  the  anode 
GDL  side  using  Eq.  (54). 


2. 2. 9. 2.  Electro-osmotic  drag.  The  protons  travelling  through 
the  pores  of  the  Nation  membrane  (from  the  anode  to  the  cath¬ 
ode)  generally  drag  one  ore  more  water  molecules  along  with 
them.  The  degree  to  which  movement  causes  water  movement 


D)  =  10  6  exp 


2416 


303 


1 

T 


+  0.0264X2  —  0.000671A3) 


^  (2.563  —  0.33A. 
for  A  >  4 


(60) 


As  for  the  drag  coefficients,  other  expression  of  the  Nation  diffu¬ 
sivity  has  been  published  [39] .  Eq.  (60)  is  valid  only  if  A  is  greater 
than  4;  if  A  is  lower  than  for  other  expressions  of  the  Nation  dif¬ 
fusivity  can  be  found  [30]  but  it  corresponds  to  working  points 
where  the  membrane  is  severely  dehydrated. 


2. 2. 9. 4.  Water  mass  balance.  Combining  Eqs.  (58)  and  (59)  the 
total  water  flux  in  Nation  is  given: 


j  M  _  n 

*'H20,  net  ~  z^drag 


7 

2F 


Pdry  _  dA 

E>x  — 
Mm  ax 


(61) 


As  proposed  in  [34,19],  the  previous  equation  can  be  rearranged 
setting  net  =  a  j/( 2  F )  where  a  is  unknown  and  denotes 
the  ratio  of  water  flux  in  the  membrane  to  hydrogen  flux: 


d  A 
dx 


Sat 


drag  22 


j 

2  F  Pdry  77^. 


(62) 
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From  Eq.  (42),  the  total  water  mass  flow  (kg  s  l)  in  Nation  is, 

^H20,net  =  MH20  ^0,net  ^tot  =  ^H20  Ot  ^tot  (63) 

In  the  ordinary  differential  equation  (62)  there  are  two 
unknowns,  X  and  a.  This  equation  cannot  be  solved  analytically 
since  DA  is  a  function  of  X.  However,  if  D,  does  not  change 
too  much  with  X  when  X  is  greater  than  4.0  (it  is  assumed  to  be 
constant  in  [30]),  an  estimation  of  D\  can  be  found  and  Eq.  (62) 
can  be  solved  analytically. 

The  water  diffusivity  in  the  Nation  is  calculated  as  follows: 


Dx  ~  D{k) 


(64) 


where  (X)  is  the  average  water  content  in  the  membrane  given 
by  Eq.  (56). 

Solving  Eq.  (62),  X  profile  in  the  membrane  can  be  found: 


A.(x)  = 


11  a 


Sat 

drag 


+  C  exp 


/  jM m«gg  \ 

y22  F  Pdry  J 


(65) 


The  two  unknowns  a  and  C  (integration  constant)  can  be  found 
from  the  two  boundary  conditions  X(0)  and  X(8U)  which  are 
computed  from  Eq.  (54). 


2. 2. 9. 5.  Membrane- specific  resistance  and  voltage  drop.  The 
Nation  conductivity  crM  (S  cm-1)  is  given  by  [34]: 

1268  ( —  -  -  |  (66) 
\303  T)\ 

where 


aM(T,  X)  =  0303  k (  A )  exp 


Or303K(A)  —  0\X  — 02 


(67) 


where  o\  —  0.005193  and  02  =  0.00326. 

The  membrane-specific  resistance  rM  (£2  cm-2)  is  obtained 
by  integrating  the  local  resistance  over  the  membrane  thickness 
8U  (cm): 


,  M 


A  dA 

')  crM(T,  X(x» 


(68) 


Finally,  after  integrating  Eq.  (68): 


Table  2 


Nafion  membrane  types  (thickness) 


Name 

<5m  (mil) 

8m  (|xm) 

Nafion  117 

7 

178 

Nafion  115 

5 

127 

Nafion  112 

2 

51 

3.  Results  and  discussion 

The  simulation  of  the  model  has  been  performed  under  both 
Smash  from  Dolphin  Integration  and  Simplorer  from  Ansoft. 
Simplorer  proposes  optimisation  functionalities  like  the  Simplex 
or  Genetic  Algorithms  which  have  been  used  to  identify  the 
model  parameters  of  the  Nexa  module  from  Ballard. 

The  identification  method  is  shown  in  Fig.  13:  experimen¬ 
tal  data  have  been  imported  in  the  simulation  tool  (2D  look-up 
tables).  Experimental  measurements  of  the  temperature  (the 
model  does  not  take  into  account  temperature  dynamics),  current 
and  air  mass  flow  feed  the  model:  experimental  and  simulated 
voltages  are  compared.  The  squared  voltage  error  is  integrated 
along  all  the  cycle:  it  is  the  objective  function  to  minimise  by 
varying  parameters  (Stou  A,  b ,  C,  V2,  VA). 

Because  the  Nexa  module  is  a  closed-loop  system  and  all  the 
data,  like  the  humidity  of  the  air,  are  not  known,  it  has  been 
assumed  that  the  cathode  and  anode  volumes  are  equal  and  that 
the  cathode  is  perfectly  humidified,  that  is,  Ph2o(Q  =  psat- 

A  parameter  sweep  analysis  has  been  performed  around  the 
solution  given  by  the  optimisation.  The  integral  quadratic  error 
has  been  plotted  vs.  each  parameter  variation  around  the  optimal 
solution  in  Fig.  14. 

As  shown  in  the  figures,  three  parameters  (activation  losses 
parameters  and  total  active  surface  area)  have  an  influence  on 
the  integral  of  the  quadratic  error.  Moreover  it  can  be  seen  that 
the  obtained  solution  is  optimal:  each  of  the  three  parameters 
obtained  by  means  of  the  optimisation  is  at  the  minimum  of  the 
error. 

However,  two  parameters  (the  concentration  losses  coeffi¬ 
cient  C  and  volumes  of  the  cathode  and  anode  V)  have  nearly  no 
influence  on  the  error.  The  explanations  are  different  for  each 
parameter: 


M 


2  exp[1268((l/r)- (1/303))] 


j  Mm( 22  0-1  a  -  2  <72  «drag) 
ln(22 eri  a  +  2o\  j Mm(-n^/22F Piry  °m) 

-  2  <72  4 ?‘g)  -  J  ft  drag  +  22  F  pdly  D{1) 


[-22  Fpdry  D 


(A) 


ln(22  <7i  a  +  2  o\  C  ftjfag  -  2  cr2  «draa2) 


MragJ 


(69) 


The  membrane  thickness  depends  on  the  membrane  type,  they 
are  summarised  in  Table  2. 

The  voltage  drop  across  the  membrane  is  simply  deduced 
from  the  membrane- specific  resistance  rM  and  the  current  den¬ 
sity  j  from  the  Ohm  law: 

^Ohm  =  ~nrMj  (70) 


•  The  concentration  losses  coefficient :  The  fuel  cell  used  to 
validate  the  model  does  not  allow  to  work  at  high  currents 
(with  possibly  low  stoichiometry  ratio)  where  the  concentra¬ 
tion  losses  become  predominant.  In  the  conditions  where  the 
tests  have  been  performed  the  concentration  losses  are  not  sig¬ 
nificant  enough  to  show  any  influence  on  the  fuel  cell  stack 
voltage.  To  identify  this  parameter  tests  should  be  performed 
at  higher  currents  and  lower  air  stream  (/.<?.,  lower  air  stoi¬ 
chiometry  ratio)  to  decrease  the  partial  pressure  of  oxygen 
inside  the  GDL.  These  tests  cannot  be  performed  on  the  Nexa 
module  because  the  system  enters  emergency  when  the  partial 
pressure  of  oxygen  drops. 

•  The  cathode  and  anode  volumes :  The  Nexa  module  works 
at  a  nearly  atmospheric  pressure  with  open-mode  cathode. 
There  is  no  pressure  valve  to  control  the  pressure  inside  the 
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Fig.  13.  Parameters  identification  method. 


cathode  gas  channels.  This  means  that  the  pressure  dynamics 
is  not  really  high  for  an  air  flow  change  and  the  pressure 
stays  roughly  around  the  atmospheric  pressure.  To  identify 
this  parameter,  tests  should  be  performed  on  a  high  pressure 
(around  2-3  bar)  fuel  cell  stack  where  air  flow  changes  have 
a  bigger  influence  on  the  pressure  dynamics. 

Fig.  15  shows  the  experimental  current  profile  which  has  been 
used  to  feed  the  model.  The  experimental  and  simulated  voltages 
are  plotted  in  Fig.  16:  it  can  be  seen  that  the  model  fits  quite  well 
to  the  experimental  results.  Some  differences  between  the  model 
and  the  simulation  are  relatively  big  at  low  current.  These  dif¬ 
ferences  could  be  explained  by  the  assumption  that  the  cathode 
channels  are  saturated:  the  Nexa  humidity  and  heat  exchanger 
time  response  is  not  known  and  can  have  an  influence  on  the 
water  content  of  the  inlet  air.  The  big  time  transients  in  the  volt¬ 
age  between  1500  s  and  2000  s  can  be  explained  by  the  time 
constant  of  the  humidifier  and/or  the  time  constant  of  the  mem¬ 
brane  water  content:  the  membrane  model  is  assumed  to  be  in 
steady  state.  A  more  refined  model  has  to  be  built  to  take  into 
account  this  phenomena. 

However,  most  of  the  voltage  errors  between  the  model  and 
experiment  appear  for  low  current:  they  do  not  have  a  big  influ¬ 
ence  on  the  simulated  power  as  seen  in  Fig.  17.  The  experimental 
and  simulated  powers  show  very  good  agreement:  this  model 


can  consequently  be  used  with  power  converters  or  other  energy 
sources  in  a  fuel  cell  system  for  power  applications  like  the  one 
from  [26]. 

The  membrane- specific  resistance  depends  on  the  tempera¬ 
ture  (Fig.  18),  the  membrane  water  content  (Fig.  19)  and  current 
density  as  seen  in  (69).  Mann  et  al.  [40]  proposed  an  empirical 
expression  for  the  specific  membrane  resistance  (see  Fig.  20)  and 
the  results  predicted  by  Eq.  (69)  show  good  agreements  with  the 
Mann’s  formula.  It  can  be  seen  in  Fig.  20  that  the  specific  resis¬ 
tance  is  not  constant  and  cannot  be  considered  constant  as  it  is 
proposed  in  a  great  deal  of  papers.  For  example,  at  the  beginning 
of  the  simulation,  the  specific  resistance  is  equal  to  80  Q  cm-2 
and  at  1500  s  where  the  current  is  the  same  (see  Fig.  15)  the 
specific  resistance  is  lower  than  60£2  cm-2  because  the  water 
content  of  the  membrane  and  the  temperatures  are  different. 

According  to  the  resistance  evolution,  the  voltage  should 
decrease  between  1700  s  and  2000  s.  However,  the  current  is 
really  low  and  the  effect  of  the  resistance  on  the  voltage  is  not 
predominant.  The  water  content  at  the  cathode  side  is  decreasing 
thanks  to  the  temperature:  the  partial  pressure  of  water  inside  the 
cathode  gas  channel  decreases  and  consequently  the  partial  pres¬ 
sure  of  oxygen  is  slightly  increasing.  According  to  the  Nernst 
equation,  if  the  oxygen  partial  pressure  increases,  the  voltage 
increases:  this  phenomenon  explains  the  small  increase  of  the 
stack  voltage  between  1700  s  and  2000  s. 
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V 

Fig.  14.  Parameters  sensitivity  analysis:  (a)  A;  (b)  b\  (c)  Sto t',  (d)  C\  (e)  V. 


Fig.  15.  Experimental  stack  current. 


Fig.  16.  Experimental  and  simulated  voltages. 
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Fig.  17.  Experimental  and  simulated  power. 


Fig.  18.  Experimental  stack  temperature. 


tent. 


Fig.  20.  Simulated  membrane-specific  resistance  with  comparison  with  resis¬ 
tance  given  by  Mann  et  al.  [40,41]. 


4.  Conclusion 

This  study  presents  a  VHDL-AMS  modelling  of  a  complex 
multi-domain  energy  conversion  system:  a  fuel  cell  stack.  A 
comparative  study  between  the  different  modelling  approaches 
(bond  graphs,  electrical  equivalent  circuits)  is  given  to  show 
the  great  advantages  of  the  VHDL-AMS  language  in  the  design 
process  of  fuel  cell  systems.  The  fuel  cell  model  includes  the 
water  transport  in  the  membrane. 

The  fuel  cell  stack  model  fits  the  experimental  results  and 
the  VHDL-AMS -based  modelling  approach  shows  its  powerful 
capabilities  on  a  fuel  cell  model,  which  is  a  complex  multi- 
domain  energy  system.  The  presented  model  is  able  to  predict 
the  voltage  and  the  power  of  the  fuel  cell  with  a  good  accuracy 
taking  into  account  the  water  content  of  the  membrane.  This 
last  point  is  really  important  to  design  the  air  supply  system 
(compressor  and  humidifier)  and  its  associated  control. 

The  modelling  approach  allows  the  design  team  to  split  the 
work  into  several  parts  (concurrent  engineering)  and  validate 
each  part  independently.  This  work  is  a  first  step  in  building  a 
complete  VHDL-AMS  model  of  a  fuel  cell  system  including 
auxiliaries  and  power  converters.  The  next  steps  will  consist 
in  improving  the  fuel  cell  stack  model  by  including  its  ther¬ 
mal  behaviour  and  the  electrodes  flooding  phenomenon,  which 
has  to  be  accounted  for  if  some  working  conditions  are  consid¬ 
ered  (high  current  density  and  low  stoichiometry  ratio).  Finally, 
the  fuel  cell  stack  auxiliaries  and  the  power  converters  will 
be  modelled.  Once  the  model  of  each  sub-part  has  been  built 
and  validated,  the  complete  system  can  be  assembled.  The  fuel 
cell  system  virtual  prototype  will  help  the  designers  determine 
the  optimal  design  of  the  system  (hybridization  with  super¬ 
capacitors,  batteries,  etc.)  and  its  optimal  control  before  its 
implementation  on  a  real  fuel  cell  system. 
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